## Joining, by = c("Eleve", "Groupe")
One row per outcome and per assessment time
# create long format data
dat_long <-
dat %>%
pivot_longer(!c(participant, experimentatrice, eleve, group, ID, age, sexe, missing, s1_bin, s1_0_1, s2_bin, s2_0_3, s3_bin, s3_0_3, s4_bin, s4_0_4), names_to = "outcome", values_to = "value") %>%
filter(!outcome %in% c("t0_corsi_env_nb", "t1_corsi_env_nb",
"t0_mdc_env_nb", "t1_mdc_env_nb"))
# identify timing of assessment
dat_long$time <- ifelse(grepl("t1", dat_long$outcome), "t1", "t0")
# remove assessment time to the name of outcomes
dat_long$outcome <- stringr::str_replace_all(dat_long$outcome, "t0_", "")
dat_long$outcome <- stringr::str_replace_all(dat_long$outcome, "t1_", "")
dat_long$outcome <- stringr::str_replace_all(dat_long$outcome, " ", "_")
dat_long <- rename_outcomes(dat_long)
dat_long$value <- ifelse(dat_long$flip == "yes",
- dat_long$value, dat_long$value)
One row per outcome (i.e., each assessment time is stored in a
different column).
We also z-scaled all the outcomes
# we translate the data + center our outcomes at t0 and t1 + assess the pre/post correlation for each of our outcomes
dat_semiwide_raw <-
dat_long %>%
pivot_wider(names_from = "time", values_from = "value") %>%
dplyr::group_by(outcome) %>%
dplyr::mutate(z0 = as.numeric(as.character(scale(t0))),
z1 = as.numeric(as.character(scale(t1))),
r = as.numeric(as.character(cor.test(~t0 + t1)$estimate))) %>%
dplyr::ungroup()
dat_semiwide <- rename_outcomes(dat_semiwide_raw)
Sum all z-scores by cognitive function / test
Demographic characteristics of our sample depending on the group.
dat_table1 <-
dat_comp %>%
group_by(group) %>%
mutate(sex = if_else(sexe == "M", 0, 1)) %>%
summarise(mean_age = mean(age),
mean_sex = mean(sex, na.rm = TRUE))
dat$sex <- factor(dat$sexe)
dependent = "group"
table1::table1(~ age + sex | group, data=dat, overall = "Overall",
render.continuous=c(.="Mean / Median<br>SD<br>[min, max]"))
| Control (N=22) |
Experimental (N=20) |
Overall (N=42) |
|
|---|---|---|---|
| age | |||
| Mean / Median SD [min, max] |
63.1 / 62.0 3.64 [58.3, 69.5] |
63.8 / 64.2 3.33 [57.8, 69.5] |
63.4 / 63.3 3.47 [57.8, 69.5] |
| sex | |||
| F | 12 (54.5%) | 11 (55.0%) | 23 (54.8%) |
| M | 10 (45.5%) | 9 (45.0%) | 19 (45.2%) |
Outcome distribution at baseline
dat_inhib_dist <- subset(dat_long, grepl("Inhibition", dat_long$GEN_outcome) & time != "t1")
ggplot(dat_inhib_dist, aes(x = value)) +
geom_histogram(aes(y=..density..), fill = "grey", alpha = 0.8) +
geom_density(alpha = 0.3, aes(fill = group)) +
facet_wrap(. ~ outcome, scales = "free") +
theme_bw() +
theme(legend.position = "bottom")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
dat_memory_dist <- subset(dat_long, grepl("memory", dat_long$GEN_outcome) & time != "t1")
ggplot(dat_memory_dist, aes(x = value)) +
geom_histogram(aes(y=..density..), fill = "grey", alpha = 0.8) +
geom_density(alpha = 0.3, aes(fill = group)) +
facet_wrap(. ~ outcome, scales = "free") +
theme_bw() +
theme(legend.position = "bottom")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
dat_EF_dist <- subset(dat_long, grepl("Executive", dat_long$GEN_outcome) & time != "t1")
ggplot(dat_EF_dist, aes(x = value)) +
geom_histogram(aes(y=..density..), fill = "grey", alpha = 0.8) +
geom_density(alpha = 0.3, aes(fill = group)) +
facet_wrap(. ~ GEN_outcome, scales = "free") +
theme_bw() +
theme(legend.position = "bottom")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
dat_othr_dist <- subset(dat_long, !grepl("memory", dat_long$GEN_outcome) & !grepl("Inhibition", dat_long$GEN_outcome) & !grepl("Executive", dat_long$GEN_outcome) & time != "t1")
ggplot(dat_othr_dist, aes(x = value)) +
geom_histogram(aes(y=..density..), fill = "grey", alpha = 0.8) +
geom_density(alpha = 0.3, aes(fill = group)) +
facet_wrap(. ~ outcome, scales = "free") +
theme_bw() +
theme(legend.position = "bottom")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
res_main = res_main %>% mutate(across(where(is.numeric), round, digits=3))
DT::datatable(res_main,
rownames = FALSE,
extensions = "Buttons",
options = list( # options
scrollX = TRUE,
dom = c('tB'),
buttons = c('copy', 'csv', 'excel', 'pdf', 'print'),
scrollY = "600px",
pageLength = 30,
order = list(3, 'asc'),
columnDefs = list(
list(width = '70px',
targets = "_all"),
list(className = 'dt-center',
targets = "_all"))))
tab.main <- data.frame(
'General outcome' = res_main$outcome,
TOST = res_main$TOST)
metaviz::viz_forest(x = res_main[, c("d", "d_se")],
variant = "classic",
col = "Greys", xlab = "SMD", annotate_CI = TRUE,
study_table = tab.main,
type = "study_only",
text_size = 2.5,
x_limit = c(-1.5, 1.5),
x_breaks = seq(-3, 3, 1))
dat_plot <-
dat_comp %>%
pivot_longer(
cols = c(z0, z1),
names_to = "time",
values_to = "value",
)
dat_plot$time[dat_plot$time == "z0"] <- "t0"
dat_plot$time[dat_plot$time == "z1"] <- "t1"
dat_inhib <- subset(dat_plot, grepl("Inhibition", dat_plot$GEN_outcome))
ggplot(dat_inhib, aes(x = time, y = value, color = group, fill = group)) +
# geom_line() +
geom_boxplot(alpha = 0.2) +
facet_grid(. ~ GEN_outcome, scales = 'free') +
theme_bw() +
theme(legend.position="bottom")
dat_memory <- subset(dat_plot, grepl("memory", dat_plot$GEN_outcome))
ggplot(dat_memory, aes(x = time, y = value, color = group, fill = group)) +
geom_boxplot(alpha = 0.2) +
facet_grid(. ~ GEN_outcome, scales = 'free') +
theme_bw() +
theme(legend.position="bottom")
dat_EF <- subset(dat_plot, grepl("Executive", dat_plot$GEN_outcome))
ggplot(dat_EF, aes(x = time, y = value, color = group, fill = group)) +
geom_boxplot(alpha = 0.2) +
facet_grid(. ~ GEN_outcome, scales = 'free') +
theme_bw() +
theme(legend.position="bottom")
dat_othr <- subset(dat_plot, !grepl("memory", dat_plot$GEN_outcome) & !grepl("Inhibition", dat_plot$GEN_outcome) & !grepl("Executive", dat_plot$GEN_outcome))
ggplot(dat_othr, aes(x = time, y = value, color = group, fill = group)) +
geom_boxplot(alpha = 0.2) +
facet_grid(. ~ GEN_outcome, scales = 'free') +
theme_bw() +
theme(legend.position="bottom")
row = length(unique(dat_semiwide$outcome))
res_S1 <- data.frame(outcome = rep(NA, row),
flip = rep(NA, row),
n_cont = rep(NA, row),
n_exp = rep(NA, row),
d = rep(NA, row) ,
d_se = rep(NA, row) ,
d_low = rep(NA, row),
d_up = rep(NA, row),
d_tot = rep(NA,row),
B = rep(NA, row),
SE = rep(NA, row),
p = rep(NA, row))
a = 0
for (out in unique(dat_semiwide$outcome)) {
# initialize some settings and re-initializing them as NA at each loop
a = a+1
mod = d = means = NA
# subset dataframe ----------------
dat_i = subset(dat_semiwide, outcome == out)
## ancova pre/post adjusted
mod = lm(z1 ~ group + z0 + age + sexe, data = dat_i)
means = emmeans::emmeans(mod, ~factor(group))
# identify critical information ------------------------------------------------------------
res_S1$outcome[a] = out
res_S1$flip[a] = unique(dat_i$flip)
res_S1$n_cont[a] = nrow(subset(dat_i, group == "Control"))
res_S1$n_exp[a] = nrow(subset(dat_i, group == "Experimental"))
res_S1$B[a] = summary(mod)$coefficients[2,1]
res_S1$SE[a] = summary(mod)$coefficients[2,2]
res_S1$p[a] = summary(mod)$coefficients[2,4]
# extract results ---------
d = smd(x = means, df = dat_i, level = 95)
res_S1$d[a] = d[[1]]
res_S1$d_low[a] = d[[2]]
res_S1$d_up[a] = d[[3]]
res_S1$d_tot[a] = d[[4]]
}
## Warning in qnorm(1 - cer): NaNs produced
## Warning in qnorm(1 - cer): NaNs produced
## Warning in qnorm(1 - cer): NaNs produced
## Warning in qnorm(1 - cer): NaNs produced
res_S1$d_se = (res_S1$d_up - res_S1$d_low) / (2*1.96)
### TOST analysis
res_S1$TOST <- NA
for (i in 1:nrow(res_S1)) {
res_S1$TOST[i] <- max(
abs(tsum_TOST(m1 = res_S1$d[i], m2 = 0, sd1 = 1, sd2 = 1,
n1 = res_S1$n_cont[i], n2 = res_S1$n_exp[i],
low_eqbound = -0.8, high_eqbound = 0.8, alpha = 0.05,
var.equal = FALSE, bias_correction = FALSE,
eqbound_type = "raw")$smd$dhigh),
abs(tsum_TOST(m1 = res_S1$d[i], m2 = 0, sd1 = 1, sd2 = 1,
n1 = res_S1$n_cont[i], n2 = res_S1$n_exp[i],
low_eqbound = -0.8, high_eqbound = 0.8, alpha = 0.05,
var.equal = FALSE, bias_correction = FALSE,
eqbound_type = "raw")$smd$dlow))
}
res_S1 <- rename_outcomes(res_S1)
res_S1 = res_S1 %>% mutate(across(where(is.numeric), round, digits=3))
DT::datatable(res_S1,
rownames = FALSE,
extensions = "Buttons",
options = list( # options
scrollX = TRUE,
dom = c('tB'),
buttons = c('copy', 'csv', 'excel', 'pdf', 'print'),
scrollY = "600px",
pageLength = 30,
order = list(3, 'asc'),
autoWidth = TRUE,
columnDefs = list(
list(width = '70px',
targets = "_all"),
list(className = 'dt-center',
targets = "_all"))))
tab.supp<- data.frame(
'General outcome' = res_S1$GEN_outcome,
Outcome = res_S1$outcome,
TOST = res_S1$TOST)
metaviz::viz_forest(x = res_S1[, c("d", "d_se")],
variant = "classic",
col = "Greys", xlab = "SMD", annotate_CI = TRUE,
group = res_S1$GEN_outcome,
study_table = tab.supp,
type = "study_only",
text_size = 2.2,
x_limit = c(-1.5, 1.5),
# N = tab.prim$N,
x_breaks = seq(-3, 3, 1))
res_S2 <- data.frame(outcome = rep(NA, 9),
flip = rep(NA, 9),
n_cont = rep(NA, 9),
n_exp = rep(NA, 9),
d = rep(NA, 9) ,
d_se = rep(NA, 9) ,
d_low = rep(NA, 9),
d_up = rep(NA, 9),
d_tot = rep(NA,9),
B = rep(NA, 9),
SE = rep(NA, 9),
p = rep(NA, 9))
a = 0
for (out in unique(dat_comp$GEN_outcome)) {
# initialize some settings and re-initializing them as NA at each loop
a = a+1
mod = d = means = NA
# subset dataframe ----------------
dat_i = subset(dat_comp, GEN_outcome == out)
## ancova pre/post adjusted
mod = lm(z1 ~ group + z0, data = dat_i)
means = emmeans::emmeans(mod, ~factor(group))
# identify critical information ------------------------------------------------------------
res_S2$outcome[a] = out
res_S2$flip[a] = unique(dat_i$flip)
res_S2$n_cont[a] = nrow(subset(dat_i, group == "Control"))
res_S2$n_exp[a] = nrow(subset(dat_i, group == "Experimental"))
res_S2$B[a] = summary(mod)$coefficients[2,1]
res_S2$SE[a] = summary(mod)$coefficients[2,2]
res_S2$p[a] = summary(mod)$coefficients[2,4]
# extract results ---------
d = smd(x = means, df = dat_i, level = 95)
res_S2$d[a] = d[[1]]
res_S2$d_low[a] = d[[2]]
res_S2$d_up[a] = d[[3]]
res_S2$d_tot[a] = d[[4]]
}
res_S2$d_se = (res_S2$d_up - res_S2$d_low) / (2*1.96)
### TOST analysis
res_S2$TOST <- NA
for (i in 1:nrow(res_S2)) {
res_S2$TOST[i] <- max(
abs(tsum_TOST(m1 = res_S2$d[i], m2 = 0, sd1 = 1, sd2 = 1,
n1 = res_S2$n_cont[i], n2 = res_S2$n_exp[i],
low_eqbound = -0.8, high_eqbound = 0.8, alpha = 0.05,
var.equal = FALSE, bias_correction = FALSE,
eqbound_type = "raw")$smd$dhigh),
abs(tsum_TOST(m1 = res_S2$d[i], m2 = 0, sd1 = 1, sd2 = 1,
n1 = res_S2$n_cont[i], n2 = res_S2$n_exp[i],
low_eqbound = -0.8, high_eqbound = 0.8, alpha = 0.05,
var.equal = FALSE, bias_correction = FALSE,
eqbound_type = "raw")$smd$dlow))
}
res_S2 = res_S2 %>% mutate(across(where(is.numeric), round, digits=3))
DT::datatable(res_S2,
rownames = FALSE,
extensions = "Buttons",
options = list( # options
scrollX = TRUE,
dom = c('tB'),
buttons = c('copy', 'csv', 'excel', 'pdf', 'print'),
scrollY = "600px",
pageLength = 30,
order = list(3, 'asc'),
columnDefs = list(
list(width = '70px',
targets = "_all"),
list(className = 'dt-center',
targets = "_all"))))
tab.S2 <- data.frame(
'General outcome' = res_S2$outcome,
TOST = res_S2$TOST)
metaviz::viz_forest(x = res_S2[, c("d", "d_se")],
variant = "classic",
col = "Greys", xlab = "SMD", annotate_CI = TRUE,
study_table = tab.S2,
type = "study_only",
text_size = 2.5,
x_limit = c(-1.5, 1.5),
x_breaks = seq(-3, 3, 1))
dat_eval = subset(dat, s4_bin == 1)
paste0(round(nrow(dat_eval)/sum(!is.na(dat$s4_bin))*100, 3), "% of the sample succeed at a beebot mastering test after the last training session, with a mean score of ", round(mean(dat_eval$s4_0_4), 3), "/4")
## [1] "89.474% of the sample succeed at a beebot mastering test after the last training session, with a mean score of 3.824/4"
dat_semiwide$z_change = dat_semiwide$z1 - dat_semiwide$z0
dat_s3 = subset(dat_semiwide, !is.na(s4_0_4))
dat_s3$s4_bin = factor(dat_s3$s4_bin)
summary(lmerTest::lmer(z_change ~ ordered(s4_0_4) + (1|participant) + (1|outcome), data = dat_s3))
## boundary (singular) fit: see ?isSingular
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: z_change ~ ordered(s4_0_4) + (1 | participant) + (1 | outcome)
## Data: dat_s3
##
## REML criterion at convergence: 967.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.1003 -0.4310 -0.0649 0.4627 3.6076
##
## Random effects:
## Groups Name Variance Std.Dev.
## participant (Intercept) 0.01782 0.1335
## outcome (Intercept) 0.00000 0.0000
## Residual 1.14216 1.0687
## Number of obs: 323, groups: participant, 19; outcome, 17
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.14038 0.11677 15.00000 1.202 0.248
## ordered(s4_0_4).L -0.15088 0.21720 15.00000 -0.695 0.498
## ordered(s4_0_4).Q -0.06353 0.23355 15.00000 -0.272 0.789
## ordered(s4_0_4).C 0.13878 0.24882 15.00000 0.558 0.585
##
## Correlation of Fixed Effects:
## (Intr) o(4_0_4).L o(4_0_4).Q
## or(4_0_4).L -0.618
## or(4_0_4).Q -0.169 -0.431
## or(4_0_4).C 0.093 -0.102 -0.398
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
res_s3_bin = lmerTest::lmer(z_change ~ s4_bin + (1|participant) + (1|outcome), data = dat_s3)
## boundary (singular) fit: see ?isSingular
emmeans(res_s3_bin, ~s4_bin)
## s4_bin emmean SE df lower.CL upper.CL
## 0 0.3562 0.1845 16.80 -0.0335 0.746
## 1 0.0254 0.0633 9.07 -0.1176 0.168
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
ggplot(dat_s3, aes(x = s4_0_4, y = z_change)) +
geom_jitter(width=0.2, size = 2, alpha = 0.4) +
theme_bw()
ggplot(dat_s3, aes(x = factor(s4_bin), y = z_change)) +
geom_jitter(width=0.1, size = 2, alpha = 0.4) +
theme_bw()
dat_agg <-
dat_long %>%
# pivot_wider(names_from = "time", values_from = "value") %>%
dplyr::group_by(outcome) %>%
dplyr::mutate(z = scale(value)) %>%
dplyr::ungroup() %>%
dplyr::group_by(time, ID, group) %>%
summarise(value_z = sum(z))
## `summarise()` has grouped output by 'time', 'ID'. You can override using the
## `.groups` argument.
ID <-
dat_agg %>%
group_by(ID) %>%
arrange(time) %>%
mutate(diff = value_z - lag(value_z, default = first(value_z))) %>%
filter(time == "t1") %>%
mutate(Improvement = if_else(diff > 0, "Improvement", "Deterioration"))%>%
dplyr::select(ID, Improvement) %>%
filter(!is.na(Improvement))
dat_line <- merge(dat_agg, ID, by = "ID")
ggplot(dat_line, aes(x = time, y = value_z, color = group, group = ID)) +
geom_point(size = 2.5, alpha = 0.4) +
geom_line(size = 1, linetype = "dotted", alpha = 0.8) +
theme_bw() +
facet_grid(group ~ Improvement) +
theme(legend.position = "none")